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^ : Abstract 

^ I We investigate both analytically and by numerical simulation the kinetics of 

^ , a microscopic model of hard rods adsorbing on a linear substrate, a model 

o, 

O , which is relevant for compaction of granular materials. The computer sim- 
ulations use an event-driven algorithm which is particularly efficient at very 

O^ I long times. For a small, but finite desorption rate, the system reaches an 

Q I equilibrium state very slowly, and the long-time kinetics display three succes- 

T-H ■ 

0^ I sive regimes: an algebraic one where the density varies as 1/t, a logarithmic 



one where the density varies as l/ln(t), followed by a terminal exponential 
approach. The characteristic relaxation time of the final regime, though incor- 
C rectly predicted by a mean field arguments, can be obtained with a systematic 

o' 

^ ' gap-distribution approach. The density fluctuations at equilibrium are also 

k><( \ investigated, and the associated time-dependent correlation function exhibits 

C^ I a power law regime followed by a final exponential decay. Finally, we show 

that denser particle packings can be obtained by varying the desorption rate 

during the process. 
68.45.Da, 61.43.-j, 64.70.Pf 
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I. INTRODUCTION 

The packing of granular materials is somewhat paradoxical. A child learns quickly that 
it is necessary to shake a bucket in order to pack the sand inside, but physicists cannot 
provide a completely satisfactory explanation of the densification process. The absence of a 
reference model, like the hard-sphere fluid for liquid-state physics or Ising model for phase 
transitions and magnetism, is at the origin of the slow progress in this field, despite a renewal 
of interest in recent years [jT| . 

To capture the main features of the packing mechanism, the experimental study of a 
model system as simple as possible can help in building a reference theory. In this spirit, 
Knight et al |^ have considered a system of monodisperse spherical beads. The packing 
process is realized by placing beads in a column which is tapped periodically with a given 
intensity. In a first series of experiments, they demonstrated that the density increases 
monotonically with the number of taps for various intensities of tapping. The very slow 
increase of density was analyzed, and a formula expressing the density in terms of the 
inverse of the logarithm of the number of taps was shown to be more accurate than any of 
the other suggestions [^. Such behavior is common to models whose geometric exclusion 
effects dictate the kinetics of densification, i.e., models in which addition of new particles is 
exponentially limited by the inverse of the free volume [§-0]. 

In a second series of experiments, Nowak et al [|],^ showed the presence of re- 
versible/irreversible cycles. The beads in an initially loosely compacted state were vibrated 
for fixed periods with a sequence of increasing vibrational intensity, causing the density to 
increase monotonically. The sequence was then reversed so that the powder was vibrated 
with decreasing intensity. The density, however, continued to increase showing that the 
initial branch is irreversible. When the initial sequence of increasing vibration was repeated, 
the second branch was retraced confirming that it is reversible. 

In the same experiments, Nowak et al [^,0 monitored the power spectrum of the density 
fluctuations around the steady state for different values of the tapping strength. A two-step 



spectrum was observed characterized by two frequencies that both increase with increasing 
tapping strength. To account for the slow kinetics of compactification, the existence of 
a "reversible" steady-state, and the fluctuation power spectrum, they proposed a simple 
adsorption-desorption or "parking lot" model. Partial analyses of this model have already 
been reported [|I^-|T^ We present here a comprehensive description of the kinetics, including 
the final exponential regime and of the fluctuations around the steady state (equilibrium). 
We first present the model in section II. We detail in section III the specific algorithm 
that we have developed for enhancing the frequency of rare events in the late stages of the 
densification process. In section IV we study the densification kinetics. By using a gap 
distribution analysis we derive an expression for the time of relaxation towards equilibrium 
and the form of the gap distribution function in the limit of small desorption; the results 
compare very well with the simulation data. A short account of this derivation has been 
given in Ref [1^. In section V, the time-dependent density-density correlation function 



is studied in the equilibrium state. The correlation function displays two well separated 
timescales, corresponding to two relaxation steps, and this can be interpreted by a simple 
model. In section VI, we show that a faster densification can be obtained by changing the 
adsorption rate during the process. 

II. THE MODEL 

In the adsorption-desorption model, particles are placed in a D-dimensional space at 
randomly selected positions with a constant rate fc+. If the trial particle does not overlap 
any previously adsorbed particle, the new particle is accepted. In addition, all adsorbed 
particles are subject to removal (desorption) at random with a constant rate k-. The 
one-dimensional version of the model, in which the substrate is a line and the objects are 
hard rods, has been solved in some limiting cases. When A;_ = 0, the adsorption is totally 
irreversible and the process corresponds to a ID Random Sequential Adsorption (RSA) for 
which the kinetics are known exactly |T^. Without a relaxation mechanism, this process is 



driven towards a non-equilibrium state and the long-time kinetics are given by an algebraic 
scaling law, p^—p{t) ~ 1/t, with p^o — 0.747.. (when the substrate is empty at the beginning 
of the process). When kj^ = 0, starting with any configuration of particles, one obtains an 
analytical solution for this uniform desorption process ||T5[. The limit fc_ -^ 0"^, which allows 
a small but non-zero possibility of rearrangement of the particles on the line, leads to a final 
density equal to 1. It is worth noting the finite discontinuity between the final density of 
this case (A;_ -^ 0^) and the RSA jamming limit (A;_ = 0). Moreover, the final density is 
independent of the initial configuration of particles on the line, whereas the jamming limit 
for the RSA process depends strongly on the initial state of the line. For A;_ — > 0+, accurate 
descriptions have been obtained |lTO| , p!T|| . In this case, the process cleanly divides into two 
sub-processes. The initial phase consists of an irreversible adsorption and it is followed 
by an infinite sequence of desorption-adsorption events in which a rod detaches from the 
surface and the gap that is created is immediately filled by one or two new rods. The latter 
possibility causes the system to evolve continuously to the close-packed state with p = I 
as |TU| , pT| 1 — p(t) ~ 1/ ln(t) where t now represents a rescaled time. For the general case, 
where both k^ and k_ are non zero, a complete solution is not available. 

The properties of the adsorption-desorption model depend only on the ratio K = k^/k^. 
A large value of K then corresponds to a small desorption rate. If time is expressed in units 
of kl^^, the densification kinetics is given by 

where $(t), the insertion probability, is the fraction of the substrate that is available for the 
insertion of a new particle. The presence of a relaxation mechanism, i.e., competing des- 
orption and adsorption with an equilibrium constant K, implies that the system eventually 
reaches a steady state that corresponds to an equilibrium configuration of hard particles with 
Peq = K^cqiPcq), where peq dcuotcs the equilibrium density. At equilibrium, the insertion 
probability is given exactly by 

$eq(p) = (l-p)exp(-p/(l-p)). (2) 



Inserting Eq. (^ in Eq. (|T]) leads to the following expression for the equilibrium density: 

P'^'^-l + LUK) ^^^ 

where L^(x) (the Lambert- W function) is the solution of x = ye"^. In the limit of small K, 
the isotherm takes the Langmuir form, p^q ~ K/{l+K), while for large K, p^g ~ 1 — 1/ ln{K). 
At small values of K, equilibrium is rapidly obtained, but at large values the densification 
process is dramatically slowed. 

III. SIMULATION ALGORITHM 

A naive method of simulating the adsorption-desorption process would attempt to ran- 
domly insert a new particle at fixed time intervals. This approach, however, is extremely 
inefficient at high densities since most attempts to add new particles are unsuccessful after 
an initial period. We have developed a general algorithm that enables us to investigate in 
detail the kinetics of the adsorption process at arbitrarily long times and for arbitrarily large 
values of the equilibrium constant, K. Instead of using a fixed time step, the algorithm is 
based on adsorption or desorption events. In this section we describe the general features 
of the algorithm that could be to simulate an adsorption-desorption process of arbitrar- 
ily shaped particles in any dimension. Later, we detail the methodology for the hard-rod 
system. 

The total rate of adsorption and desorption events is 

RUt) = m+pit)/K. (4) 

The quiesence, or waiting time, is the time interval between any two successive events that 
alter the state of the system. Let F{t) denote the probability that the waiting time is 
greater than r. Since successive events are considered to be independent, 

F{t + At) = F(r)(l - i?totAr + OiAr^)). (5) 

Taking the limit At ^ we obtain 



F(r)=exp(-i?totr). (6) 

A uniformly distributed random number, < .^i < 1, may be used to sample a random 
waiting time consistent with this distribution: 

r = -ln(ei)/i?tot. (7) 

Once the quiesence time has been selected, the nature of the event is determined stochasti- 
cally by defining 

and choosing a second uniformly distributed random number, < ^2 < 1- If C2 < f'dit) 
the event at time t is a desorption and a randomly selected particle is removed from the 
system. If ^2 > i^d{t) the event is an adsorption and a new particle is placed randomly in the 
available surface. We have validated the method by applying it to the Langmuir equation 
where $ = 1 — p/pmax and the kinetics, as well as the isotherm, are known exactly. This 
method is quite general and can apply to a range of adsorption-desorption processes. 

In the simulation, the initial state of the system is an interval of length L (measured in rod 
lengths) bounded by two immovable rods centered at positions xq = —1, xn+i = L+1. For an 
arbitrary configuration of A^ additional rods, whose centers are located at {xj, i = 1, . . . N}, 
the total available length is known exactly: Lq = J2f=^^ max(xi+i — Xj — 2, 0). At each step 
of the simulation, the total rate of adsorption and desorption events is determined from 
R = Lq + N/K. A waiting time is sampled from the exponential waiting time distribution 
using Eq.(|^ and the type of event is decided with Eq. (H). If the event is adsorption, a new 
particle is placed in the available length. The probability that a particular gap is occupied is 
equal to its available length divided by the total available lengthLo. Thus a random number 
^3 is generated and the position of the particle on the available length is ^3-^0, which means 
that the gap between particles j and j + 1 is occupied where j is defined by the following 
equation: 



i-1 j 

Y^ max(xi+i - Xi - 2, 0) < ^sLq < Yl max(xj+i - Xj - 2, 0). (9) 

j=0 i=0 

Note that the adsorption event is uniform and is always accepted. If the event is desorption 
a particle is selected at random and removed from the surface. Note that the desorption 
probability is independent of the length of time that the particle has been on the sur- 
face. The available line is updated; it always decreases (increases) following an adsorption 
(desorption) event. The simulation procedure thus generates a sequence of configurations 
(ti, Ni, Lqi), (^2, ^^2, -^02)5 ••• and one knows the state of the system at an arbitrary time, t. 
To insure good statistics, several thousands of independent simulations must be run for each 
value of the desorption rate. We used system lengths, L, from 400 to 5000. 

IV. GAP DENSITY APPROACH 

The adsorption-desorption model can be described in terms of gap distribution functions. 
The one-gap distribution function, G{h,t), represents the density of voids of length h; the 
time evolution of G{h, t) is given by 

^^i^ = -H(h - l)(h - l)G(h,t) + 2 r dh'G(h',t) 
ot Jh+l 

-^G{h, t) + ^^^ jy dh'Gih', h-l-h',t), (10) 

where G{h, h', t) is the two-gap distribution function associated with the probability of find- 
ing two neighboring gaps (separated by one particle) of length h and h' , and H{x) is the 
Heaviside step function {H{x) = 1 for x > l,H{x) = otherwise). The first two terms on 
the right-hand side of Eq. (p!0D correspond to loss and gain terms due to adsorption while 
the remaining two are due to desorption. Similarly, the time evolution of G{h, h', t) is given 
by 

^^^^'^^ = -{H{h -l){h-l) + H{h' - l)ih' - l))Gih, h', t) 

/"CO /"OO 

+ / dh"G{h",h',t)+ dh"G{h,h",t)+G{h + h' + l,t) 

Jh+l Jh'+l 



_ ±G{h, h', t) + ^ , '' / dh"G{h\ h-l- h", h', t) 
K Kp[t) Jo 

+ ^y~^^ r '' dh"G{h, h'-i- h", h\ t), (11) 

Kp[t) Jo 
where G{h, h', h", t) is the three-gap distribution function. The kinetics of the process is thus 
given by an infinite hierarchy of equations involving an infinite set of muhi-gap distribution 
functions. 

The quantities of interest can be expressed in terms of integrals of the one-gap distribu- 
tion function. In particular, the insertion probability $(t) is given by 



oo 



<l>(t) = / dh{h-l)G{h,t), (12) 

and we have the following sum rules: 

/•oo 

p{t) = / dhG{h,t), (13) 

Jo 

and 

/•oo 

l-p(t) = / dhhG{h,t). (14) 

Jo 

(One also has G{h, t) = /q°° dh'G{h, h', t) = /q°° dh'G{h', h, t) and similar integrals for higher- 
order terms). The steady-state solution of Eq. ([To| ) is known and corresponds to the equi- 
librium hard-rod system with a gap distribution function given by |T0|,11 



and all higher-order distribution functions satisfying the factorization property, 

Gcq{hi, /l2, ■ ■ ■ , hn, p) = Geq{hi, p)Geq(/i2, p) " " " Gcq{hn, p) ■ (16) 

In order to find a solution for the kinetics of the process, one must truncate the hierarchy 
by means of a closure ansatz. The simplest closure is provided by by an adiabatic (mean- 
field) treatment. There, one assumes that, at any density p(t), the structure of the adsorbate, 
as characterized by the gap distribution functions, is that of an equilibrium system at density 
Pit). 
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V. DENSIFICATION KINETICS (AT CONSTANT K » I) 

A. A succession of regimes 

We focus here on the small desorption limit {K >> 1). For an initially empty line, there 
are three different successive kinetic regimes. The first stage is dominated by adsorption 
events and the process displays an RSA-like behavior, characterized by a 1/t power law 
dependence. For p{t) > 0.7, adsorption becomes slower and desorption can no longer be 
ignored, which allows particle rearrangements on the line and, eventually, insertion of addi- 
tional particles. The densification mechanism requires the rearrangement of an increasing 
number of particles in order to open a hole large enough for the insertion an additional parti- 
cle. The process is similar to what occurs in the limit k ^ 0^, and the kinetics is dominated 
by a l/ln(t) behavior |Ty,|TT|. For large but finite values of K, this densification regime 
continues until the density is very close to the equilibrium (steady-state) value, pe^^K). In 
the final regime, the desorption term becomes comparable to the adsorption term, and an 
exponential approach to equilibrium is observed. Figure ^ illustrates the three successive 
regimes. 

B. Exponential approach to equilibrium 

The exponential regime is illustrated in Figure ^. In an adiabatic (mean-field) treatment 
the insertion probability, $(t), satisfies an equation similar to Eq. (^ with p{t) in place of 
Peg- Denoting 6p{t) = p{t) — p^o, with p^o = Peq{K), one obtains from Eq. (|lD, at first order 
in density, 

j^6p = -Tmf{K)6p + 0{6p^) (17) 



with 



Imf[K) = — , (18j 

~ ln{Ky/K when K is large, (19) 



which is equivalent to a relaxation time given by K/\n.{KY for large K. 

In Fig. ^, the relaxation rate is plotted as a function of K: the dashed curve gives 
the mean- field prediction, Eq. (|19[), and open circles correspond to the best exponential fit 
to the simulation results. It is evident that the mean-field analysis gives a poor estimate 
of the relaxation rate for large K. This failure can be understood by the noting that the 
mean-field assumption leads to a characteristic time for the rearrangement of $ of the order 
if/ln(i^)^, i.e., much smaller than K, the characteristic time for desorption. Since in the 
absence of surface diffusion process, significant rearrangement can only occur on a timescale 
longer than K the system is unable to adjust rapidly enough in order to change significantly 
the available surface function on a timescale of order i^/ln(i^)^. In Fig 0,we display the 
insertion probability, $(p), for several large values of K] it is worth noting that $ first 
follows the RSA curve until it reaches a value close to the equilibrium one at which point 
it plateaus and evolves very weakly towards equilibrium. The process clearly deviates from 
the adiabatic description in which the insertion probably is given, at all densities, by the 
equilibrium curve, <l>eq(p)- 

We now turn to a description in terms of the gap distribution approach. To obtain the 
leading term in the exponential approach towards equilibrium, when K is very large (but 
finite), we assume that, as for the steady state (or equilibrium), \G{h,t)\ ~ exp(— II/i), with 
n ~ Ini^ ~ (1 — p)^^ (see Eqs. (|15|) and (§), when K is very large). As a consequence, if 
one defines p„(t) = /„"+^ G{h, t)dh and $„(t) = /"+^(/i - l)G(/i, t)dh, then p„ ~ $„ ~ iT"", 
so that if one looks for the dominant behavior in 1/ K, it is sufficient to consider the first 
intervals in h. As in the adiabatic approximation, one can expand the gap densities in power 
of 5p{t) and keep only the linear term which gives rise to the exponential decay. In the final 
regime, where the density is close to the steady state, we first assume that the deviation of 
the gap distribution function from its equilibrium form, 6G{h,t) = G{h,t) — Geq{h), can be 
expressed as an expansion in Sp{t) where only the first term is kept. Let us then denote 

d\nG{h,p) 



A{h) = p^ ^^ 
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(20) 

Poo 



and 



A{h, h') = pc 



d In G{h,h',p) 



dp 



(21) 



Poo 



After inserting Eq. (^) in Eq. ([To|) A{h) can be rewritten for < /i < 1 as 



2-7 



A(/i) 









(22) 



where Pqo = Poo/(1 — Poo) is the dimensionless equihbrium pressure for poo = Peg and 



7 = irr = -i^(5p(t)/5p(t)| 



Poo- 



(23) 



From Eqs. (|l|) and (|T2D, one obtains 

7 = 1 - p2 / dhhe-^^^A{l + /i). 

JO 

whereas the sum rules in Eq. (|l^) and (|l^ give, respectively, 

("OO 

Poo / dhe-^^^Aih) = 1 
Jo 

- Poo / dhhe-^^^Aih) = 1. 
Jo 



(24) 



(25) 



(26) 



When integrating the two sides of Eqs. ( p2D between and 1, one obtains 

(2 - 7) = 2P^ / dhhe-^^^A{l + h) + 0{1/K). (27) 

^ u 

Combining Eq. ( P^ ) with Eq. ( |27D yields 7 = 0(l/i^). Thus, the relaxation rate F goes 
essentially as l/K"^ instead of the 1/K dominant behavior predicted by the mean-field treat- 
ment. In order to have a more explicit expression for 7 and F, it is necessary to calculate 
the integral on the right-hand side of Eq. (l2^) to O^l/K"^): 



(2-7) 



1-^f/ e-^-^A(l + /i) + A(l) 
K \Jo 



( P3 fl 

2-2 0-- ■ 



l\h + l)e-^^^ A{2 + h)\ 
+ 0(1/7^2), (28) 



which leads to 
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^K = 2Poo (a(0) + PooA(l) - Pi J\h + 1)6-''-'' A{2 + h)^ + 0{1/K). (29) 

An explicit expression for •yK, Eq. ( ^9]) thus requires the knowledge of the gap distribution 
function for 3 > h > 1. The kinetic equation for the gap distribution function when h > 1 
is then rewritten by inserting Eqs. (pOD and (^) in Eq. ([lOD , which gives 

(h-l + ^^) A{h) = -{h-l)+ f dh'A{h', h-l-h') 

+ £^ / d/ie-^-e^ -"M(l + h). (30) 

Combining Eq. ( P^ ) and Eq. ( |50D for /i = 1, one finally gets 

7/r = 2Poo (-4(0) - P^ I' rf/i/ie-^°°'^A(2 + /i)) + 0{l/K). (31) 

Since the system evolves close to equilibrium, we further assume that the factorization 
property for the two-gap distribution function is valid to 0{1/K), i.e., A{h,h') = A{h) + 
A{h') + 0{1/K) [|TB[. Eqs. (PDj ) and (P^) then become a closed set of equations for A{h) 
to 0(1/7^). The solution is given in the Appendix, as well as the explicit expression for 
T{K) = •y/K, Eq. ( |A17| ). As an illustration, the leading terms of T{K) in powers of In(i^) 



are obtained as 

(InKf (InKV (InK) ^ 1 , 
r - 2^^ - 4^^ + 2^^ + 0(— ). (32) 

The prediction of Eq. ( |A17| ), shown as the full curve in Fig. |^, gives a good agreement with 
the results obtained from an exponential fit to the simulation data, whereas the mean-field 
predictions fail completely (dashed curve) (the dotted curve corresponds to the first term of 
the right hand- side of Eq. (^)). 

VI. FLUCTUATIONS AROUND EQUILIBRIUM 

For times much larger than the relaxation time, the density no longer evolves (on aver- 
age), but fluctuates around its equilibrium value. Note that in this regime, the fluctuation- 
dissipation theorem and the time translational invariance are both valid. We have calculated 
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the time-dependent correlation function C(t) of the density fluctuations, Sp(t) = pit) — poo, 
around the equilibrium state. Starting from the 2— time correlation function, 

C(i' + M') = <K« + «W)>-<P('' + ')><K«')>, (33) 

^ ' ^ < pit'Y > - < pit') >2 ' ^ ^ 

we have numerically verified that when t' » l/F, C(t' + t, t') becomes time translationnally 
invariant, i.e. C(t + t',t') = C(t). (Conversely, when 1 << t' « l/F, one observes aging 
phenomena [0, but we postpone the discussion of this phenomena to a future publication). 
Because of the very long relaxation time, we found that the calculation of the correlation 
function is more efficient by using Eq. ( p3D instead of taking the usual time average on 
a single system ||18[ . Results from the simulation are shown in Fig.^ for two large values 



of K. At short and intermediate times, the decay of C{t) is non-exponential, whereas 
at long times the kinetics follows an exponential decay. The two regimes, or relaxation 
steps, can be interpreted as follows: the first consists of a "fast" adsorption-desorption 
process without appreciable densification of the system, whereas the second corresponds to 
the linear-response regime and, as predicted by Onsager's regression hypothesis, it shows 
the same final exponential dependence as the final approach of p{t) towards poo in the 
densification process. In this second relaxation step, C{t) ~ e^'"* where F is given by 
Eqs. (H) and ([MTI) . 

Close to equilibrium and for large values of K, the adsorption and desorption events 
can be considered as spatially uncorrected, and the system can be represented as a set of 
independent gaps in which a particle is adsorbed or not. This assumption does not account 
for rearrangements which occur at long times, but is valid for short times. When a particle 
is adsorbed, the gap is characterized by the distribution geq{h) which is the distribution 
probability of finding a particle such the total length of right and left gap is equal to h, i.e. 

9e,ih) = — f dh'G,,ih')G,,ih - h') 

Poo -^0 

= p^PlheM-Pooli). (34) 

Once the particle has desorbed, the gap is characterized by the distribution Geqih + 1). The 
two distributions are calculated at equilibrium and their ratio is given exactly by: 

13 



9e,{h) Kh- ^ ^ 

For a given gap of length {h + 1), we have a two-state (particle-hole) stochastic process in 
which the rates, associated to desorption and adsorption respectively, are 1/K and h. The 
average probability for having a particle in the gap is equal to 

For a given gap of length {h + 1), the (unnormalized) correlation function of the density 
fluctuations Ch(t) due to the two-state stochastic process is given by |^ 



Cnit) = -(^^^fM-{yK + h)t). (37) 

Assuming that the adsorption-desorption events giving rise to the two-state process only 
seldom affect simultaneously two neighboring particles, one can write the (unnormalized) 
correlation function as a superposition of correlation functions occurring in parallel in the 
different gaps, weighted by the distribution geq{h), i.e.. 



oo 



Cshort = / dhge,{h)Ch{t) (38) 







At equilibrium, the variance of the density fluctuations can be calculated exactly 

<(M'>=Poo(l-Poo)', (40) 

so that, the normalized correlation function C{t) at short times, C short, can be written as 

which reduces to a power law, 1/t, when In(i^) << t << K. This result is equivalent to the 
ln(u;)-behavior already predicted along similar lines similar by Kolan et al [jl9|. The insets 
in Figs. ^,b illustrate the excellent agreement between Eq. (^) and the simulation data. 

It is worth noting that the 1/t behavior is reminiscent of the pure RSA asymptotic regime, 
where it occurs as a consequence of the filling of small isolated pieces of the available fraction 
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of the line, whose lengths go to zero when t -^ oo. With similar arguments, one thus expects 
to have in higher dimensions a t'"^/^ behavior, leading to an cj'^^^^/'^) power-law dependence 
for the power spectrum. In particular, this predicts a power law tu"^/^ for D = 2, which is 
compatible with the experimental data in vibrated granular media ||^. In dimensions higher 
than 1, our prediction differs from that of Kolan et al [|r^ since their analysis leads to a 
1/uj dependence in the power spectrum. A numerical study of the two-dimensional version 
of the adsorption-desorption model should settle this point. 

VII. DENSIFICATION REGIME AND MULTISTEP PROCESS 

The very slow-exponential approach to equilibrium with T{K) = 0(1/ K'^) when K is 
very large implies that $(p) increases with p when p is sufficiently large (see section VB and 
the Appendix). Since in the first (RSA-like) regime $(p) decreases, there always exists a 
density pm where (9<l>(p)/9p|p„^ = 0. Fig. || displays a log-log plot of $ as a function of time 
for various values of K. One notice (i) that pm is an increasing function of K and (ii) that 
the minimum of $ is always very close but smaller than the equilibrium value, ^eql-^), smd 
smaller than p/K, which is due to the fact that the density is an increasing function of time. 

In Fig. § the density is plotted as a function of time for different values of K. The 
curves on the left part of the figure correspond to an adiabatic process where the available 
surface function is replaced by the equilibrium formula, Eq. (H), and which corresponds to a 
process where rapid diffusion on the substrate allows a instantaneous equilibration after each 
desorption and adsorption event. For all values of K, the adiabatic process is much faster 
than the corresponding adsorption-desorption model process. Moreover, for an adiabatic 
process, the density is at all times always a monotonically increasing function of K. For the 
adsorption-desorption model, on the other hand, the density is not always monotonic in K. 
In Fig 1^, for example, the system with K = 500 has a higher density than the system with 
K = 1000 for 4 ^ ln(t) ^ 8. The existence of a minimum in $ is a sufficient condition for 
this phenomenon. It follows that for a given finite time, the densification can be made more 
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efficient by changing the desorption rate during the process. Fig. |^ compares the densities 
obtained by using either a single value oi K = 1000 or a sequence of varying K, starting 
from 1000 at t = 0, passing through a minimum, and finishing at the same value K = 1000 
when t = 1000 . One clearly observes that a larger final density is reached with the multistep 
process. Such a phenomenon, which is also the source of the reversible-irreversible cycles 
observed by Nowaket al [^, has been already observed and quantitatively analyzed in an 
irreversible adsorption model [|l5l. However, the determination of the optimum densification 



strategy, which has significant applications to vibratory compaction of granular materials, 
is still an open problem. 

ACKNOWLEDGEMENTS 

J.T. thanks the National Science Foundation for financial support. The Laboratoire 
de Physique Theorique des Liquides is Unite Mixte de Recherche No 7600 au CNRS. The 
Laboratoire de Physique Theorique is Unite Mixte de Recherche No 8627 au CNRS. 

APPENDIX: SOLUTION FOR THE GAP DISTRIBUTION FUNCTION IN THE 

LINEAR RESPONSE REGION 

For convenience, we introduce the notation B{h) = A{h) — 1, B{h,h') = A{h,h') — 2, 
etc. For < /i < 1, by taking into account that 7 = 0{1/K), Eq. (UT) can be reexpressed 



as 



5(/i)e-^-^ - 5(1)% = P^ /' dh'e~^-'^'B{l + h') + 0{e-^-^^ / K) (Al) 

K JO 



and for larger gaps, Eq. (^) can be rewritten for /i > 0, as 
2 fh 2P 



(h + —)B(1 + h)= dh'B(h', h-h') + — l^e^-'^ / U^dh'e-'^^'' B(2 + h') + 0(1/K^). 
K Jo K Jh 

(A2) 



Assuming that the factorization property is valid to a 0{1/K), i.e., 
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B{h, h') = B{h) + B{h') + ^%^ + 0(1/7^2)^ 

K 



(A3) 



with < /i < 1, and 0{1/K) or 0{1/K'^) designate functions that are uniformly of order 
1/ K or 1/ K'^ on the interval [0, 1], we can derive from Eq. (^^) 

2 rh — l 1 rh roo 

(h + —)B(l + h)-2 dh'Bih') = — / dh'C(h',h-h') + 2P^ dh' e-^^^^ -^^ Bi2 + h' 

K Jo K [jo Jh 

+ 0{1/K'), (A4) 



and 



r-l + h 



(A5) 



(l + h)B(2 + h)-2 f dh'B(h')=0(l/K). 
Jo 

When /I >> 2/i^, Eq. (p^ simplifies to 

/ifi(l + h)-2 f dh'B{h') = 0{l/K). (A6) 

io 

one gets the 



(A7) 



.Pcxj/l' 



Deriving Eq. ( [A^ ) with respect to /i and inserting the result in Eq. 
following differential equation, 



/5- J^(^^(^ + ^)) - ih^^^^^ + ^^^ + ^^^^^^ + ^^^ = ^^^/^^ 



with l>h» 2/ K^ whose solution is 



5(1 + h) = b{2 - P^h) + c 



P^{h + 2/K) 



+ (2 - Poo/i)e-^°° / 

JO 



/i' + 2/K 



(A8) 



where 6 and c are constants and 1 > /i >> 2/K. The corresponding solution for B{h) is 
then 



S(^) = 6(1 - Pooh) + c 



-Poo{l-h) 



+ (1 - Poo/i)e-^°° /' 

JO 



dh'- 



■,Pac:h' 



(A9) 



/o h' + 2/ir_ 

It is important to stress that the above equations give the solution only when h>2/K. To 
satisfy Eqs. ( [Al| ) and (|A3| ) when h ~ 2/K or smaller, on must include in B{1 + h) an 
additional component that is a 0(1/^) when h » 2/K and is non negligible only when 
h ~ 2/K. The full solution for < /i < 1 is then obtained as 
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5(1 + h) = 6(2 - P^h) + c 
d 



"^ ^ + (2 - P^/.)e-^^ 



+ 



P^{h + 2/ir) 
0(l/ir). 



dh'- 



.Pooh' 



/i' + 2/K_ 

(AlO) 



fs:(/i + 2/ir) 

It is easy to verify that Eq. ( |X9| ) is still the full solution to a 0(l/i^) for < h < 1 and 
that, from Eq. ([A5|) , the solution for 5(2 + /i), < /i < 1, is equal to : 

1 



5(2 + /i) 



^^ ^ ^^ ^5(2) + bh{4 - P^h) -^p- Pooh)e-''-^'-'^ 



+ (2-4Poo/i+(Poo/i)2)e^^-|^ - 



h' gPoah' 



+ Oil/K) 



(All) 



+ 2/K_ 

The constants 6, c, d, and -8(2) are determined by the various sum rules as well as by the 
condition, which comes from the structure of the hierarchy of kinetic equations, that B{h) 
is a piecewise continuous function, namely. 



B{1) = 2b-^ + ^ = b{l-P^) + c[l + (1 - Poo)e-^-E,(Poo)] + 0{1/K), 
B{2) = 6(2 - Poo) + c[--^ + 1(2 - Poo)e-^-E,(Poo)] + 0{l/K), 

where Ei{x) = e^/x + e"" J^ dtexp{-t)/{x - t)^. 
The result can be expressed as 



(A12) 
(A13) 



P(0) 



P^ + 1 



P(1) = (P«, + 1) 

P(2) = (Poo + 1) 



-(Poo-l)+2(Poo + l 

2 - Poo + 2 



;i-Poo$-)(i + (i-Poo)$,^ 



1 + (1 - Poo)($f + $r) - Poo{2 - Poo)$f $ 
(Poo + 1) (1-Poo$r)(^oo-1+Poo(2-Poo)$ 



r) 



Poo 1 + (1 - Poo)($f + $r) - Poo{2 - Poo)$f $° 

where we have introduced the notation $°° = e~^°°Ei{Pao) and $f^ = e^°°-Ei(Poo) with 
Ei{x) = J^ dt^-j—. The values of c and d can be trivially derived from the above equations. 



(AM) 
(A15) 

(A16) 



The relaxation rate T = ^K can be obtained by inserting the above solution into Eq. (^ 
which leads to 



ri^2_7ir = 2Poo(i + P« 



1 + t/o 



(A17) 



with 






When Pf^ ^ oo, one has 



^r-^(l-;^ + ^ + 0(A)], (A19) 

-'oo\ -foo -'oo -'00/ 

^r-;^(l + ;^ + ;^ + 0(-^))- (A20) 

-Too \ -Too -Too -"^oo / 



Inserting Eqs. ( [XTP| ) and (^) in Eqs.(^) and (|M7|) leads to Eq. (| 
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FIGURES 
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FIG. 1. Linear-logarithmic plot of the adsorbed density as a function of time for a large value 
of K (K=5000). The process is characterized by three slow kinetic regimes: (i) RSA-like regime 
whose final stage is described by a 1/t behavior, (ii) l/ln(i) regime, and (iii) exponential approach 
towards equilibrium. 
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FIG. 2. (a) Final exponential approach of the density p to its equilibrium value peq for two 
large values of K. (b) Relaxation rate for the approach to equilibrium F versus K. Upper curve: 
prediction from mean-field approximation, Eq. (^). Dotted curve: leading tem of Eq. (32). Full 
curve, Eq. ( A17] ) . Open circles: best exponential fit to the numerical simulations. 
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FIG. 3. The insertion probability $ as a function of density p for various large values of K 
{K = 100,500, 1000,5000). The dashed curve corresponds to a process without desorption (RSA 
process) and the dotted curve corresponds to the equilibrium insertion probability, Eq. (P). 
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FIG. 4. Equilibrium density-density correlation function C{t) versus time for (a) K = 500 and 
(b) K = 1000. The inset in the upper right corner displays the first step in the decay of the corre- 
lation function (full curve) as well as the predicted short-time formula, Cshor'ti't)^ Eq. (^) (dashed 
curve). The other inset shows the exponential decay of C{t) at long times on a logarithmic-linear 
plot. 
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FIG. 5. Log-log plot of the insertion probability <1> as a function of time for various large values 
of K {K = 100,500, 1000,5000). Note that for K > 100, ^ displays a minimum which is smaller 
than the equilibrium value, (dotted lines) 
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FIG. 6. Linear-Log plot of the density versus time for different values of K. The left curves 
correspond to the adiabatic process. The right curves correspond to the adsorption-desorption 
model. Notice that in the latter case the curves for different values of K always cross, a phenomenon 
absent in the adiabatic process 
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FIG. 7. Density increase over 1500 time units for a process with a single value of K {K = 1000) 
and for multistep process in which the sequence of K is shown on a linear-logarithmic plot in the 
inset (After t = 1000, K stays constant). Note the large enhancement of packing efficiency in 
the latter process. This effect is absent in the adiabatic approximation where the density is a 
monotonically increasing function of K. 
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